Q1: Plot surface temperature Climatology for
import numpy as np
import matplotlib.pyplot as plt
import xarray as xr
import cartopy as cp
import cartopy.feature as cfeature
import cartopy.crs as ccrs
from pylab import *
may = xr.open_dataset('../monsoon/datasets/skt.may.nc')
# may = may.skt
clim_may = may.mean(dim='time')
from cartopy.util import add_cyclic_point
clim_may_data = clim_may['skt']
lon = clim_may.coords['lon']
print("Original shape -", clim_may_data.shape)
lon_idx = clim_may_data.dims.index('lon')
wrap_clim_may, wrap_lon = add_cyclic_point(clim_may_data.values, coord=lon, axis=lon_idx)
print("New shape -", wrap_clim_may.shape)
# Plotting the data
# lon = clim_may.lon
lat = clim_may.lat
fig = plt.figure(figsize=(16,16))
ax = fig.add_subplot(1,1,1,projection=ccrs.PlateCarree(central_longitude=0.0, globe=None))
mp = ax.imshow(wrap_clim_may-273.5,extent=(wrap_lon.min(),wrap_lon.max(),lat.min(),lat.max()),cmap='jet',origin='upper')
plt.title('Long Term Mean Temperature (May)')
# plt.legend(['Temp'])
states_provinces = cfeature.NaturalEarthFeature(
category='cultural',
name='admin_1_states_provinces_lines',
scale='10m',
facecolor='none')
# ax.add_feature(cfeature.BORDERS,edgecolor='blue')
# ax.add_feature(states_provinces, edgecolor='blue')
ax.add_feature(cfeature.LAND)
ax.add_feature(cfeature.COASTLINE)
ax.add_feature(cfeature.OCEAN)
cbar = fig.colorbar(mp, shrink=0.3,label='Temperature (degC)')
cbar.minorticks_on()
#adding the long lat grids and enabling the tick labels
gl = ax.gridlines(draw_labels=True,alpha=0.1)
gl.top_labels = False
gl.right_labels = False
Original shape - (94, 192) New shape - (94, 193)
from cartopy.util import add_cyclic_point
months = ['may','june','july','august','september','jjas']
# lon = clim_may.lon
lat = clim_may.lat
# plotting all the clim_data graphs at once
for i in months:
clim_data = xr.open_dataset('../monsoon/datasets/skt.'+i+'.nc')
clim_data = clim_data.mean(dim='time')
# adding the cyclic points to remove the missing values in the center of the plot
data = clim_data['skt']
lon = clim_data.coords['lon']
print("Original shape -", data.shape)
lon_idx = data.dims.index('lon')
clim_wrap_data, wrap_lon = add_cyclic_point(data.values, coord=lon, axis=lon_idx)
print("New shape -", clim_wrap_data.shape)
# clim_data = clim_data.skt
# clim_data = clim_data.mean(dim='time')
fig = plt.figure(figsize=(16,16))
ax = fig.add_subplot(1,1,1,projection=ccrs.PlateCarree(central_longitude=0.0, globe=None))
mp = ax.imshow(clim_wrap_data-273.5,extent=(wrap_lon.min(),wrap_lon.max(),lat.min(),lat.max()),cmap='jet',origin='upper')
plt.title('Mean Surface Temperature ('+i.capitalize()+')')
ax.add_feature(cfeature.LAND)
ax.add_feature(cfeature.COASTLINE)
ax.add_feature(cfeature.OCEAN)
cbar = fig.colorbar(mp, shrink=0.3,label='Temperature (degC)')
cbar.minorticks_on()
gl = ax.gridlines(draw_labels=True,alpha=0.1)
gl.top_labels = False
gl.right_labels = False
plt.savefig('../monsoon/figures/Surface Temperature ('+i+').png')
Original shape - (94, 192) New shape - (94, 193) Original shape - (94, 192) New shape - (94, 193) Original shape - (94, 192) New shape - (94, 193) Original shape - (94, 192) New shape - (94, 193) Original shape - (94, 192) New shape - (94, 193) Original shape - (94, 192) New shape - (94, 193)
# function to add cyclic points
def add_cyclic(data, lon):
"""Add a cyclic point to a data array.
Parameters
----------
data : array-like
An n-dimensional array of data to add a cyclic point to.
lon : array-like
An array which specifies the coordinate values for
`data`'s longitude coordinate.
Returns
-------
The data array with a cyclic point added.
"""
cyclic_data, cyclic_lon = add_cyclic_point(data, coord=lon, axis=-1)
return cyclic_data, cyclic_lon
# Doing the same operation as above for MSLP dataset
months = ['may','june','july','august','september','jjas']
# lon = clim_may.lon
lat = clim_may.lat
# plotting all the clim_data graphs at once
for i in months:
clim_data = xr.open_dataset('../monsoon/datasets/mslp.'+i+'.nc')
# adding the cyclic points to remove the missing values in the center of the plot using the above defined function
clim_data = clim_data.mean(dim='time')
cyclic_data, cyclic_lon = add_cyclic_point(clim_data['mslp'].values, coord=clim_data.coords['lon'], axis=-1)
# clim_data = clim_data.mslp
# clim_data = clim_data.mean(dim='time')
fig = plt.figure(figsize=(16,16))
ax = fig.add_subplot(1,1,1,projection=ccrs.PlateCarree(central_longitude=0.0, globe=None))
mp = ax.imshow(cyclic_data/100,extent=(cyclic_lon.min(),cyclic_lon.max(),lat.min(),lat.max()),cmap='jet',origin='upper')
plt.title('Long Term Mean MSLP ('+i.capitalize()+')')
ax.add_feature(cfeature.LAND)
ax.add_feature(cfeature.COASTLINE)
ax.add_feature(cfeature.OCEAN)
cbar = fig.colorbar(mp, shrink=0.3,label='MSLP (hPa)')
cbar.minorticks_on()
gl = ax.gridlines(draw_labels=True,alpha=0.1)
gl.top_labels = False
gl.right_labels = False
plt.savefig('../monsoon/figures/MSLP ('+i+').png')
# Plotting the time series of skin temp and MSLP for a particular location
clim_skt = xr.open_dataset('../monsoon/datasets/skt.ymonmean.nc')
clim_mslp = xr.open_dataset('../monsoon/datasets/mslp.ymonmean.nc')
clim_mslp = clim_mslp.mslp
clim_skt = clim_skt.skt
# Plotting the time series of skin temp and MSLP for a particular location
clim_skt_land = clim_skt.sel(lat=20,lon=80,method='nearest')
clim_mslp_land = clim_mslp.sel(lat=20,lon=80,method='nearest')
clim_skt_ocean = clim_skt.sel(lat=-20,lon=80,method='nearest')
clim_mslp_ocean = clim_mslp.sel(lat=-20,lon=80,method='nearest')
# plotting double y-axis plot with mslp on one side and skin temp on the other
fig, ax1 = plt.subplots(figsize=(16,8))
ax2 = ax1.twinx()
ax1.plot(clim_skt_land-273.5,'b.--',label='Skin Temp')
ax2.plot(clim_mslp_land/100,'r.--',label='MSLP')
ax1.set_xlabel('Time')
ax1.set_ylabel('Skin Temp (degC)')
ax2.set_ylabel('MSLP (hPa)')
ax1.legend(loc='upper left')
ax2.legend(loc='upper right')
plt.title('Time Series of Skin Temp and MSLP at Land (20N,80E)')
plt.savefig('../monsoon/figures/Time Series of Skin Temp and MSLP at Land (20N,80E).png')
# adding a trendline to both the curves and plotting fitting function alongwith the difference of mslp and skin temperature
# fitting a linear function to the data
import numpy as np
from scipy.interpolate import make_interp_spline
# Plotting the time series of skin temp and MSLP for a particular location
clim_skt = xr.open_dataset('../monsoon/datasets/skt.ymonmean.nc')
clim_mslp = xr.open_dataset('../monsoon/datasets/mslp.ymonmean.nc')
clim_mslp = clim_mslp.mslp
clim_skt = clim_skt.skt
# Plotting the time series of skin temp and MSLP for a particular location
clim_skt_land = clim_skt.sel(lat=20,lon=80,method='nearest')
clim_mslp_land = clim_mslp.sel(lat=20,lon=80,method='nearest')
clim_skt_ocean = clim_skt.sel(lat=-20,lon=80,method='nearest')
clim_mslp_ocean = clim_mslp.sel(lat=-20,lon=80,method='nearest')
# create a smoother set of x values
xnew = np.linspace(0, len(clim_skt_land), 300)
# interpolate the skin temperature and MSLP data using a spline function
skt_spline = make_interp_spline(range(len(clim_skt_land)), clim_skt_land-273.5)
mslp_spline = make_interp_spline(range(len(clim_mslp_land)), clim_mslp_land/100)
# calculate the interpolated values using the new set of x values
skt_smooth = skt_spline(xnew)
mslp_smooth = mslp_spline(xnew)
# plot the smoothed curves
fig, ax1 = plt.subplots(figsize=(16,8))
ax2 = ax1.twinx()
ax1.plot(xnew, skt_smooth, 'b-', label='Skin Temp')
ax2.plot(xnew, mslp_smooth, 'r-', label='MSLP')
ax1.set_xlabel('Time')
ax1.set_ylabel('Skin Temp (degC)')
ax2.set_ylabel('MSLP (hPa)')
ax1.legend(loc='upper left')
ax2.legend(loc='upper right')
plt.title('Time Series of Skin Temp and MSLP at Land (20N,80E)')
plt.savefig('../monsoon/figures/Time Series of Skin Temp and MSLP at Land (20N,80E).png')
# Importing required libraries
import numpy as np
import pandas as pd
import xarray as xr
import matplotlib.pyplot as plt
from scipy.interpolate import make_interp_spline, BSpline
# Plotting the time series of skin temp and MSLP for a particular location
clim_skt = xr.open_dataset('../monsoon/datasets/skt.ymonmean.nc')
clim_mslp = xr.open_dataset('../monsoon/datasets/mslp.ymonmean.nc')
clim_mslp = clim_mslp.mslp
clim_skt = clim_skt.skt
# Plotting the time series of skin temp and MSLP for a particular location
clim_skt_land = clim_skt.sel(lat=20,lon=80,method='nearest')
clim_mslp_land = clim_mslp.sel(lat=20,lon=80,method='nearest')
clim_skt_ocean = clim_skt.sel(lat=-20,lon=80,method='nearest')
clim_mslp_ocean = clim_mslp.sel(lat=-20,lon=80,method='nearest')
# plotting double y-axis plot with mslp on one side and skin temp on the other
fig, ax1 = plt.subplots(figsize=(16,8))
ax2 = ax1.twinx()
ax1.plot(clim_skt_land-273.5,'b.--',label='Skin Temp')
ax2.plot(clim_mslp_land/100,'r.--',label='MSLP')
# # plotting the difference of mslp and skin temp
# plt.plot(clim_skt_land+273.5-clim_mslp_land/100 ,'g.--',label='MSLP - Skin Temp')
ax1.set_xlabel('Time')
ax1.set_ylabel('Skin Temp (degC)')
ax2.set_ylabel('MSLP (hPa)')
ax1.legend(loc='upper left')
ax2.legend(loc='upper right')
# Adding smooth curves to the plot
# Define x and y values for skin temp
x1 = np.arange(len(clim_skt_land.time))
y1 = clim_skt_land.values - 273.5
# Define x and y values for MSLP
x2 = np.arange(len(clim_mslp_land.time))
y2 = clim_mslp_land.values / 100
# Smooth curves for skin temp and MSLP
s1 = make_interp_spline(x1, y1, k=3)
s2 = make_interp_spline(x2, y2, k=3)
# Plot the smooth curves
x1_new = np.linspace(x1.min(), x1.max(), 300)
x2_new = np.linspace(x2.min(), x2.max(), 300)
ax1.plot(x1_new, s1(x1_new), 'b-', label='Skin Temp smooth curve')
ax2.plot(x2_new, s2(x2_new), 'r-', label='MSLP smooth curve')
# Add actual data points to the plot
ax1.scatter(x1, y1, s=10, alpha=0.5, color='b', label='Skin Temp data points')
ax2.scatter(x2, y2, s=10, alpha=0.5, color='r', label='MSLP data points')
ax1.legend(loc='center left', frameon=False)
ax2.legend(loc='center right', frameon=False)
plt.title('Time Series of Skin Temp and MSLP at Land (20N,80E)')
plt.savefig('../monsoon/figures/Time Series of Skin Temp and MSLP at Land (20N,80E).png')
plt.show()
# Importing required libraries
import numpy as np
import pandas as pd
import xarray as xr
import matplotlib.pyplot as plt
from scipy.interpolate import make_interp_spline, BSpline
# Plotting the time series of skin temp and MSLP for a particular location
clim_skt = xr.open_dataset('../monsoon/datasets/skt.ymonmean.nc')
clim_mslp = xr.open_dataset('../monsoon/datasets/mslp.ymonmean.nc')
clim_mslp = clim_mslp.mslp
clim_skt = clim_skt.skt
# Plotting the time series of skin temp and MSLP for a particular location
clim_skt_land = clim_skt.sel(lat=20,lon=80,method='nearest')
clim_mslp_land = clim_mslp.sel(lat=20,lon=80,method='nearest')
clim_skt_ocean = clim_skt.sel(lat=-20,lon=80,method='nearest')
clim_mslp_ocean = clim_mslp.sel(lat=-20,lon=80,method='nearest')
# plotting double y-axis plot with mslp on one side and skin temp on the other
fig, ax1 = plt.subplots(figsize=(16,8))
ax2 = ax1.twinx()
ax1.plot(clim_skt_ocean-273.5,'b.--',label='Surface Temperature')
ax2.plot(clim_mslp_ocean/100,'r.--',label='MSLP')
ax1.set_xlabel('Time')
ax1.set_ylabel('Skin Temp (degC)')
ax2.set_ylabel('MSLP (hPa)')
ax1.legend(loc='upper left')
ax2.legend(loc='upper right')
# Adding smooth curves to the plot
# Define x and y values for skin temp
x1 = np.arange(len(clim_skt_ocean.time))
y1 = clim_skt_ocean.values - 273.5
# Define x and y values for MSLP
x2 = np.arange(len(clim_mslp_ocean.time))
y2 = clim_mslp_ocean.values / 100
# Smooth curves for skin temp and MSLP
s1 = make_interp_spline(x1, y1, k=3)
s2 = make_interp_spline(x2, y2, k=3)
# Plot the smooth curves
x1_new = np.linspace(x1.min(), x1.max(), 300)
x2_new = np.linspace(x2.min(), x2.max(), 300)
ax1.plot(x1_new, s1(x1_new), 'b-', label='Skin Temp smooth curve')
ax2.plot(x2_new, s2(x2_new), 'r-', label='MSLP smooth curve')
# Add actual data points to the plot
ax1.scatter(x1, y1, s=10, alpha=0.5, color='b', label='Skin Temp data points')
ax2.scatter(x2, y2, s=10, alpha=0.5, color='r', label='MSLP data points')
ax1.legend(loc='center left', frameon=True)
ax2.legend(loc='upper right', frameon=True)
plt.title('Time Series of Skin Temp and MSLP over Ocean (20S,80E)')
plt.savefig('../monsoon/figures/Time Series of Skin Temp and MSLP over Ocean (20S,80E).png')
plt.show()
Air Temperature at
for may, june and july
month = ['may', 'june', 'july']
level = ['1000','900','850']
# looping through each month and level and plotting the data for each month and level
for i in month:
for j in level:
clim_month = xr.open_dataset('../monsoon/datasets/air.'+i+'.nc')
clim_month = clim_month.mean(dim='time')
clim_month = clim_month.sel(level=j,method='nearest')
clim_data, lon = add_cyclic_point(clim_month.air, coord=clim_month.lon)
lat = clim_month.lat
fig = plt.figure(figsize=(16,16))
ax = fig.add_subplot(1,1,1,projection=ccrs.PlateCarree(central_longitude=0.0, globe=None))
mp = ax.imshow(clim_data-273,extent=(lon.min(),lon.max(),lat.min(),lat.max()),cmap='jet',origin='upper')
plt.title('Mean air temperature ('+i.capitalize()+')' + ' ' + j + 'hPa')
ax.add_feature(cfeature.LAND)
ax.add_feature(cfeature.COASTLINE)
ax.add_feature(cfeature.OCEAN)
cbar = fig.colorbar(mp, shrink=0.3,label='Temperature (degC)')
cbar.minorticks_on()
gl = ax.gridlines(draw_labels=True,alpha=0.1)
gl.top_labels = False
gl.right_labels = False
plt.savefig('../monsoon/figures/Mean air temperature ('+i+')'+ ' ' +j +'hPa.png')
plt.show()
Trying and seeing whether adding the cyclic points via interpolation works or not ??
from cartopy.util import add_cyclic_point
may = xr.open_dataset('../monsoon/datasets/skt.may.nc')
# may = may.skt
clim_may = may.mean(dim='time')
## adding the cyclic points to remove the missing data at the center of the plot
data = clim_may['skt']
lon = clim_may.coords['lon']
print("Original shape -", data.shape)
lon_idx = data.dims.index('lon')
wrap_data, wrap_lon = add_cyclic_point(data.values, coord=lon, axis=lon_idx)
print("New shape -", wrap_data.shape)
# Plotting the data
# lon = clim_may.lon
lat = clim_may.lat
fig = plt.figure(figsize=(16,16))
ax = fig.add_subplot(1,1,1,projection=ccrs.PlateCarree(central_longitude=0.0, globe=None))
mp = ax.imshow(wrap_data-273.5,extent=(wrap_lon.min(),wrap_lon.max(),lat.min(),lat.max()),cmap='jet',origin='upper')
plt.title('Long Term Mean Temperature (May)')
# plt.legend(['Temp'])
states_provinces = cfeature.NaturalEarthFeature(
category='cultural',
name='admin_1_states_provinces_lines',
scale='10m',
facecolor='none')
# ax.add_feature(cfeature.BORDERS,edgecolor='blue')
# ax.add_feature(states_provinces, edgecolor='blue')
ax.add_feature(cfeature.LAND)
ax.add_feature(cfeature.COASTLINE)
ax.add_feature(cfeature.OCEAN)
cbar = fig.colorbar(mp, shrink=0.3,label='Temperature (degC)')
cbar.minorticks_on()
#adding the long lat grids and enabling the tick labels
gl = ax.gridlines(draw_labels=True,alpha=0.1)
gl.top_labels = False
gl.right_labels = False
Original shape - (94, 192) New shape - (94, 193)
Plot low level jet stream (at 850 hPa pressure level) using NCEP U and V wind data for
Climatology for the period of 1980-2020
uwnd_may = xr.open_dataset('/home/shiv/Documents/GitHub/EES405/datasets/uwnd.may.850.nc')
uwnd_may = uwnd_may.uwnd
uwnd_may
<xarray.DataArray 'uwnd' (time: 1, level: 1, lat: 73, lon: 144)>
[10512 values with dtype=float32]
Coordinates:
* time (time) datetime64[ns] 2000-05-01
* lon (lon) float32 0.0 2.5 5.0 7.5 10.0 ... 350.0 352.5 355.0 357.5
* lat (lat) float32 90.0 87.5 85.0 82.5 80.0 ... -82.5 -85.0 -87.5 -90.0
* level (level) float32 850.0
Attributes: (12/13)
standard_name: eastward_wind
long_name: Monthly U-wind on Pressure Levels
units: m/s
cell_methods: time: mean (monthly from 6-hourly values)
precision: 2
GRIB_id: 33
... ...
var_desc: u-wind
dataset: NCEP/DOE AMIP-II Reanalysis (Reanalysis-2) Monthly Averages
level_desc: Pressure Levels
statistic: Individual Obs
parent_stat: Other
actual_range: [-58.910004 103.159996]uwnd_may = uwnd_may.squeeze(dim=['time','level'])
uwnd_may
<xarray.DataArray 'uwnd' (lat: 73, lon: 144)>
[10512 values with dtype=float32]
Coordinates:
time datetime64[ns] 2000-05-01
* lon (lon) float32 0.0 2.5 5.0 7.5 10.0 ... 350.0 352.5 355.0 357.5
* lat (lat) float32 90.0 87.5 85.0 82.5 80.0 ... -82.5 -85.0 -87.5 -90.0
level float32 850.0
Attributes: (12/13)
standard_name: eastward_wind
long_name: Monthly U-wind on Pressure Levels
units: m/s
cell_methods: time: mean (monthly from 6-hourly values)
precision: 2
GRIB_id: 33
... ...
var_desc: u-wind
dataset: NCEP/DOE AMIP-II Reanalysis (Reanalysis-2) Monthly Averages
level_desc: Pressure Levels
statistic: Individual Obs
parent_stat: Other
actual_range: [-58.910004 103.159996]uwnd_may.shape
(73, 144)
from cartopy.util import add_cyclic_point
# may = xr.open_dataset('../monsoon/datasets/skt.may.nc')
# may = may.skt
# clim_may = may.mean(dim='time')
## adding the cyclic points to remove the missing data at the center of the plot
# data = clim_may['uwnd']
lon = uwnd_may.coords['lon']
print("Original shape -", uwnd_may.shape)
lon_idx = data.dims.index('lon')
wrap_data_uwnd_may, wrap_lon = add_cyclic_point(uwnd_may.values, coord=lon, axis=lon_idx)
print("New shape -", wrap_data.shape)
# Plotting the data
# lon = uwnd_may.lon
lat = uwnd_may.lat
fig = plt.figure(figsize=(16,16))
ax = fig.add_subplot(1,1,1,projection=ccrs.PlateCarree(central_longitude=0.0, globe=None))
mp = ax.imshow(wrap_data_uwnd_may,extent=(wrap_lon.min(),wrap_lon.max(),lat.min(),lat.max()),cmap='jet',origin='upper')
plt.title('Long Term Mean Temperature (May)')
# plt.legend(['Temp'])
states_provinces = cfeature.NaturalEarthFeature(
category='cultural',
name='admin_1_states_provinces_lines',
scale='10m',
facecolor='none')
# ax.add_feature(cfeature.BORDERS,edgecolor='blue')
# ax.add_feature(states_provinces, edgecolor='blue')
ax.add_feature(cfeature.LAND)
ax.add_feature(cfeature.COASTLINE)
ax.add_feature(cfeature.OCEAN)
cbar = fig.colorbar(mp, shrink=0.3,label='Wind Speed')
cbar.minorticks_on()
#adding the long lat grids and enabling the tick labels
gl = ax.gridlines(draw_labels=True,alpha=0.1)
gl.top_labels = False
gl.right_labels = False
Original shape - (73, 144) New shape - (94, 193)
vwnd_may = xr.open_dataset('/home/shiv/Documents/GitHub/EES405/datasets/vwnd.may.850.nc')
vwnd_may = vwnd_may.vwnd
vwnd_may = vwnd_may.squeeze(dim=['time','level'])
from cartopy.util import add_cyclic_point
# may = xr.open_dataset('../monsoon/datasets/skt.may.nc')
# may = may.skt
# clim_may = may.mean(dim='time')
## adding the cyclic points to remove the missing data at the center of the plot
# data = clim_may['uwnd']
lon = uwnd_may.coords['lon']
print("Original shape -", uwnd_may.shape)
# Add cyclic point to uwnd and vwnd
lon_idx = uwnd_may.dims.index('lon')
wrap_data_uwnd_may, wrap_lon = add_cyclic_point(uwnd_may.values, coord=lon, axis=lon_idx)
wrap_data_vwnd_may,wrap_lon = add_cyclic_point(vwnd_may.values, coord=lon, axis=lon_idx)
# Compute vector addition
wind_may = np.sqrt(wrap_data_uwnd_may**2 + wrap_data_vwnd_may**2)
wind_may_direction = np.arctan2(wrap_data_vwnd_may, wrap_data_uwnd_may)
u_wind_add = wind_may * np.cos(wind_may_direction)
v_wind_add = wind_may * np.sin(wind_may_direction)
print("New shape -", wrap_data_uwnd_may.shape)
# Plotting the data
# lon = uwnd_may.lon
lat = uwnd_may.lat
fig = plt.figure(figsize=(16,16))
ax = fig.add_subplot(1,1,1,projection=ccrs.PlateCarree(central_longitude=0.0, globe=None))
# Create 2D grid of longitudes and latitudes
lon_2d, lat_2d = np.meshgrid(wrap_lon, lat)
# Compute vector addition of u-wind and v-wind
u_wind_add = u_wind_add + v_wind_add
v_wind_add = v_wind_add + u_wind_add
# Compute wind speed from u-wind and v-wind
wind_may = np.sqrt(u_wind_add**2 + v_wind_add**2)
# Plot vector addition
skip = (slice(None, None, 1), slice(None, None, 4))
ax.quiver(lon_2d[skip], lat_2d[skip], u_wind_add[skip], v_wind_add[skip], wind_may[skip],
transform=ccrs.PlateCarree(), cmap='coolwarm', scale=300)
# Plot wind speed as contour plot
mp = ax.imshow(wind_may, extent=(wrap_lon.min(),wrap_lon.max(),lat.min(),lat.max()), cmap='jet', origin='upper')
states_provinces = cfeature.NaturalEarthFeature(
category='cultural',
name='admin_1_states_provinces_lines',
scale='10m',
facecolor='none')
# ax.add_feature(cfeature.BORDERS,edgecolor='blue')
# ax.add_feature(states_provinces, edgecolor='blue')
ax.add_feature(cfeature.LAND)
ax.add_feature(cfeature.COASTLINE)
ax.add_feature(cfeature.OCEAN)
cbar = fig.colorbar(mp, shrink=0.3,label='Wind Speed')
cbar.minorticks_on()
#adding the long lat grids and enabling the tick labels
gl = ax.gridlines(draw_labels=True,alpha=0.1)
gl.top_labels = False
gl.right_labels = False
Original shape - (73, 144) New shape - (73, 145)
# making all the plots at once for may, june, july, august, september and jjas
months = ['may','june','july','august','september','jjas']
for i in months:
uwnd = xr.open_dataset('/home/shiv/Documents/GitHub/EES405/datasets/uwnd.'+i+'.850.nc')
vwnd = xr.open_dataset('/home/shiv/Documents/GitHub/EES405/datasets/vwnd.'+i+'.850.nc')
uwnd = uwnd.uwnd
vwnd = vwnd.vwnd
uwnd = uwnd.squeeze(dim=['time','level'])
vwnd = vwnd.squeeze(dim=['time','level'])
from cartopy.util import add_cyclic_point
lon = uwnd_may.coords['lon']
print("Original shape -", uwnd_may.shape)
# Add cyclic point to uwnd and vwnd
lon_idx = uwnd_may.dims.index('lon')
wrap_data_uwnd_may, wrap_lon = add_cyclic_point(uwnd_may.values, coord=lon, axis=lon_idx)
wrap_data_vwnd_may,wrap_lon = add_cyclic_point(vwnd_may.values, coord=lon, axis=lon_idx)
# Compute vector addition
wind_may = np.sqrt(wrap_data_uwnd_may**2 + wrap_data_vwnd_may**2)
wind_may_direction = np.arctan2(wrap_data_vwnd_may, wrap_data_uwnd_may)
u_wind_add = wind_may * np.cos(wind_may_direction)
v_wind_add = wind_may * np.sin(wind_may_direction)
print("New shape -", wrap_data_uwnd_may.shape)
# Plotting the data
# lon = uwnd_may.lon
lat = uwnd_may.lat
fig = plt.figure(figsize=(16,16))
ax = fig.add_subplot(1,1,1,projection=ccrs.PlateCarree(central_longitude=0.0, globe=None))
# Create 2D grid of longitudes and latitudes
lon_2d, lat_2d = np.meshgrid(wrap_lon, lat)
# Compute vector addition of u-wind and v-wind
u_wind_add = u_wind_add + v_wind_add
v_wind_add = v_wind_add + u_wind_add
# Compute wind speed from u-wind and v-wind
wind_may = np.sqrt(u_wind_add**2 + v_wind_add**2)
# Plot vector addition
skip = (slice(None, None, 1), slice(None, None, 4))
ax.quiver(lon_2d[skip], lat_2d[skip], u_wind_add[skip], v_wind_add[skip], wind_may[skip],
transform=ccrs.PlateCarree(), cmap='coolwarm', scale=300)
# Plot wind speed as contour plot
mp = ax.imshow(wind_may, extent=(wrap_lon.min(),wrap_lon.max(),lat.min(),lat.max()), cmap='jet', origin='upper')
plt.title('Wind Speed and Direction for '+i + ' at 850 hPa')
states_provinces = cfeature.NaturalEarthFeature(
category='cultural',
name='admin_1_states_provinces_lines',
scale='10m',
facecolor='none')
# ax.add_feature(cfeature.BORDERS,edgecolor='blue')
# ax.add_feature(states_provinces, edgecolor='blue')
ax.add_feature(cfeature.LAND)
ax.add_feature(cfeature.COASTLINE)
ax.add_feature(cfeature.OCEAN)
cbar = fig.colorbar(mp, shrink=0.3,label='Wind Speed')
cbar.minorticks_on()
#adding the long lat grids and enabling the tick labels
gl = ax.gridlines(draw_labels=True,alpha=0.1)
gl.top_labels = False
gl.right_labels = False
plt.savefig('/home/shiv/Documents/GitHub/EES405/monsoon/figures/'+i+'850hPa.png',dpi=300)
Original shape - (73, 144) New shape - (73, 145) Original shape - (73, 144) New shape - (73, 145) Original shape - (73, 144) New shape - (73, 145) Original shape - (73, 144) New shape - (73, 145) Original shape - (73, 144) New shape - (73, 145) Original shape - (73, 144) New shape - (73, 145)
uwnd_data = xr.open_dataset('/home/shiv/Documents/GitHub/EES405/datasets/uwnd.may.nc')
vwnd_data = xr.open_dataset('/home/shiv/Documents/GitHub/EES405/datasets/vwnd.may.nc')
uwnd_data = uwnd_data.uwnd
vwnd_data = vwnd_data.vwnd
uwnd_data = uwnd_data.sel(level=200)
vwnd_data = vwnd_data.sel(level=200)
uwnd_data = uwnd_data.squeeze(dim=['time'])
vwnd_data = vwnd_data.squeeze(dim=['time'])
vwnd_data = vwnd_data.sel(lat=slice(0,90))
<xarray.DataArray 'vwnd' (lat: 73, lon: 144)>
[10512 values with dtype=float32]
Coordinates:
time datetime64[ns] 2000-05-01
* lon (lon) float32 0.0 2.5 5.0 7.5 10.0 ... 350.0 352.5 355.0 357.5
* lat (lat) float32 90.0 87.5 85.0 82.5 80.0 ... -82.5 -85.0 -87.5 -90.0
level float32 200.0
Attributes: (12/13)
standard_name: northward_wind
long_name: Monthly V-wind on Pressure Levels
units: m/s
cell_methods: time: mean (monthly from 6-hourly values) Monthly Averages
precision: 2
GRIB_id: 34
... ...
var_desc: v-wind
dataset: NCEP/DOE AMIP-II Reanalysis (Reanalysis-2) Monthly Averages
level_desc: Pressure Levels
statistic: Individual Obs
parent_stat: Other
actual_range: [-69.100006 69.21 ]uwnd_data = xr.open_dataset('/home/shiv/Documents/GitHub/EES405/datasets/uwnd.may.nc')
vwnd_data = xr.open_dataset('/home/shiv/Documents/GitHub/EES405/datasets/vwnd.may.nc')
uwnd_data = uwnd_data.uwnd
vwnd_data = vwnd_data.vwnd
uwnd_data = uwnd_data.sel(level=200)
vwnd_data = vwnd_data.sel(level=200)
uwnd_data = uwnd_data.squeeze(dim=['time'])
vwnd_data = vwnd_data.squeeze(dim=['time'])
lon = uwnd_data.coords['lon']
print("Original shape -", uwnd_data.shape)
# Add cyclic point to uwnd and vwnd
lon_idx = uwnd_data.dims.index('lon')
wrap_data_uwnd_data, wrap_lon = add_cyclic_point(uwnd_data.values, coord=lon, axis=lon_idx)
wrap_data_vwnd_data, wrap_lon = add_cyclic_point(vwnd_data.values, coord=lon, axis=lon_idx)
# Compute vector addition
wind_may = np.sqrt(wrap_data_uwnd_data**2 + wrap_data_vwnd_data**2)
wind_may_direction = np.arctan2(wrap_data_vwnd_data, wrap_data_uwnd_data)
u_wind_add = wind_may * np.cos(wind_may_direction)
v_wind_add = wind_may * np.sin(wind_may_direction)
print("New shape -", wrap_data_uwnd_may.shape)
# Plotting the data
# lon = uwnd_may.lon
lat = uwnd_may.lat
fig = plt.figure(figsize=(16,16))
ax = fig.add_subplot(1,1,1,projection=ccrs.PlateCarree(central_longitude=0.0, globe=None))
# Create 2D grid of longitudes and latitudes
lon_2d, lat_2d = np.meshgrid(wrap_lon, lat)
# Compute vector addition of u-wind and v-wind
u_wind_add = u_wind_add + v_wind_add
v_wind_add = v_wind_add + u_wind_add
# Compute wind speed from u-wind and v-wind
wind_may = np.sqrt(u_wind_add**2 + v_wind_add**2)
# Plot vector addition
skip = (slice(None, None, 3), slice(None, None, 4))
ax.quiver(lon_2d[skip], lat_2d[skip], u_wind_add[skip], v_wind_add[skip], wind_may[skip],
transform=ccrs.PlateCarree(), cmap='coolwarm', scale=500)
# Plot wind speed as contour plot
mp = ax.imshow(wind_may, extent=(wrap_lon.min(),wrap_lon.max(),lat.min(),lat.max()), cmap='jet', origin='upper')
states_provinces = cfeature.NaturalEarthFeature(
category='cultural',
name='admin_1_states_provinces_lines',
scale='10m',
facecolor='none')
# ax.add_feature(cfeature.BORDERS,edgecolor='blue')
# ax.add_feature(states_provinces, edgecolor='blue')
ax.add_feature(cfeature.LAND)
ax.add_feature(cfeature.COASTLINE)
ax.add_feature(cfeature.OCEAN)
cbar = fig.colorbar(mp, shrink=0.3,label='Wind Speed')
cbar.minorticks_on()
#adding the long lat grids and enabling the tick labels
gl = ax.gridlines(draw_labels=True,alpha=0.1)
gl.top_labels = False
gl.right_labels = False
Original shape - (73, 144) New shape - (73, 145)
uwnd_data = xr.open_dataset('/home/shiv/Documents/GitHub/EES405/datasets/uwnd.may.nc')
vwnd_data = xr.open_dataset('/home/shiv/Documents/GitHub/EES405/datasets/vwnd.may.nc')
uwnd_data = uwnd_data.uwnd
vwnd_data = vwnd_data.vwnd
uwnd_data = uwnd_data.sel(level=200)
vwnd_data = vwnd_data.sel(level=200)
uwnd_data = uwnd_data.squeeze(dim=['time'])
vwnd_data = vwnd_data.squeeze(dim=['time'])
lon = uwnd_data.coords['lon']
print("Original shape -", uwnd_data.shape)
# Add cyclic point to uwnd and vwnd
lon_idx = uwnd_data.dims.index('lon')
wrap_data_uwnd_data, wrap_lon = add_cyclic_point(uwnd_data.values, coord=lon, axis=lon_idx)
wrap_data_vwnd_data, wrap_lon = add_cyclic_point(vwnd_data.values, coord=lon, axis=lon_idx)
# Compute vector addition
wind_may = np.sqrt(wrap_data_uwnd_data**2 + wrap_data_vwnd_data**2)
wind_may_direction = np.arctan2(wrap_data_vwnd_data, wrap_data_uwnd_data)
u_wind_add = wind_may * np.cos(wind_may_direction)
v_wind_add = wind_may * np.sin(wind_may_direction)
print("New shape -", wrap_data_uwnd_may.shape)
# Plotting the data
# lon = uwnd_may.lon
lat = uwnd_may.lat
fig = plt.figure(figsize=(16,16))
ax = fig.add_subplot(1,1,1,projection=ccrs.PlateCarree(central_longitude=0.0, globe=None))
# Create 2D grid of longitudes and latitudes
lon_2d, lat_2d = np.meshgrid(wrap_lon, lat)
# Compute vector addition of u-wind and v-wind
u_wind_add = u_wind_add + v_wind_add
v_wind_add = v_wind_add + u_wind_add
# Compute wind speed from u-wind and v-wind
wind_may = np.sqrt(u_wind_add**2 + v_wind_add**2)
# Find the indices of the maximum wind speed
max_idx = np.unravel_index(np.argmax(wind_may), wind_may.shape)
# Plot a red dot at the location of maximum wind speed
ax.plot(lon_2d[max_idx], lat_2d[max_idx], 'ko', markersize=10, transform=ccrs.PlateCarree())
# Plot vector addition
skip = (slice(None, None, 3), slice(None, None, 4))
q = ax.quiver(lon_2d[skip], lat_2d[skip], u_wind_add[skip], v_wind_add[skip], wind_may[skip],
transform=ccrs.PlateCarree(), cmap='coolwarm', scale=500)
# Plot wind speed as contour plot
mp = ax.imshow(wind_may, extent=(wrap_lon.min(),wrap_lon.max(),lat.min(),lat.max()), cmap='jet', origin='upper')
# Add colorbar
cbar = fig.colorbar(mp, shrink=0.3,label='Wind Speed')
cbar.minorticks_on()
# Add title
ax.set_title('May Wind Speed at 200 hPa')
#adding the long lat grids and enabling the tick labels
gl = ax.gridlines(draw_labels=True,alpha=0.1)
gl.top_labels = False
gl.right_labels = False
# Add land, coastlines, and ocean features
ax.add_feature(cfeature.LAND)
ax.add_feature(cfeature.COASTLINE)
ax.add_feature(cfeature.OCEAN)
plt.show()
Original shape - (73, 144) New shape - (73, 145)
# plot all the months together in one code
months = ['may','june','july','august','september','jjas']
for i in months:
uwnd_data = xr.open_dataset('/home/shiv/Documents/GitHub/EES405/datasets/uwnd.'+i+'.nc')
vwnd_data = xr.open_dataset('/home/shiv/Documents/GitHub/EES405/datasets/vwnd.'+i+'.nc')
uwnd_data = uwnd_data.uwnd
vwnd_data = vwnd_data.vwnd
uwnd_data = uwnd_data.sel(level=200)
vwnd_data = vwnd_data.sel(level=200)
uwnd_data = uwnd_data.squeeze(dim=['time'])
vwnd_data = vwnd_data.squeeze(dim=['time'])
lon = uwnd_data.coords['lon']
print("Original shape -", uwnd_data.shape)
# Add cyclic point to uwnd and vwnd
lon_idx = uwnd_data.dims.index('lon')
wrap_data_uwnd_data, wrap_lon = add_cyclic_point(uwnd_data.values, coord=lon, axis=lon_idx)
wrap_data_vwnd_data, wrap_lon = add_cyclic_point(vwnd_data.values, coord=lon, axis=lon_idx)
# Compute vector addition
wind_may = np.sqrt(wrap_data_uwnd_data**2 + wrap_data_vwnd_data**2)
wind_may_direction = np.arctan2(wrap_data_vwnd_data, wrap_data_uwnd_data)
u_wind_add = wind_may * np.cos(wind_may_direction)
v_wind_add = wind_may * np.sin(wind_may_direction)
print("New shape -", wrap_data_uwnd_may.shape)
# Plotting the data
# lon = uwnd_may.lon
lat = uwnd_may.lat
fig = plt.figure(figsize=(16,16))
ax = fig.add_subplot(1,1,1,projection=ccrs.PlateCarree(central_longitude=0.0, globe=None))
# Create 2D grid of longitudes and latitudes
lon_2d, lat_2d = np.meshgrid(wrap_lon, lat)
# Compute vector addition of u-wind and v-wind
u_wind_add = u_wind_add + v_wind_add
v_wind_add = v_wind_add + u_wind_add
# Compute vector addition of u-wind and v-wind
u_wind_add = u_wind_add + v_wind_add
v_wind_add = v_wind_add + u_wind_add
# Compute wind speed from u-wind and v-wind
wind_may = np.sqrt(u_wind_add**2 + v_wind_add**2)
# Find the indices of the maximum wind speed
max_idx = np.unravel_index(np.argmax(wind_may), wind_may.shape)
# Plot a red dot at the location of maximum wind speed
ax.plot(lon_2d[max_idx], lat_2d[max_idx], 'ko', markersize=10, transform=ccrs.PlateCarree())
# Plot vector addition
skip = (slice(None, None, 3), slice(None, None, 4))
q = ax.quiver(lon_2d[skip], lat_2d[skip], u_wind_add[skip], v_wind_add[skip], wind_may[skip],
transform=ccrs.PlateCarree(), cmap='coolwarm', scale=500)
# Plot wind speed as contour plot
mp = ax.imshow(wind_may, extent=(wrap_lon.min(),wrap_lon.max(),lat.min(),lat.max()), cmap='jet', origin='upper')
# Add colorbar
cbar = fig.colorbar(mp, shrink=0.3,label='Wind Speed')
cbar.minorticks_on()
# Add title
ax.set_title(i +' Wind Speed at 200 hPa')
#adding the long lat grids and enabling the tick labels
gl = ax.gridlines(draw_labels=True,alpha=0.1)
gl.top_labels = False
gl.right_labels = False
# Add land, coastlines, and ocean features
ax.add_feature(cfeature.LAND)
ax.add_feature(cfeature.COASTLINE)
ax.add_feature(cfeature.OCEAN)
plt.show()
Original shape - (73, 144) New shape - (73, 145)
Original shape - (73, 144) New shape - (73, 145)
Original shape - (73, 144) New shape - (73, 145)
Original shape - (73, 144) New shape - (73, 145)
Original shape - (73, 144) New shape - (73, 145)
Original shape - (73, 144) New shape - (73, 145)